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LONG-TERM  GOALS 

The  goals  of  this  work  are  to  obtain  better  understanding  of  sediment  mobilization,  transport,  and 
deposition  across  the  wave  bottom-boundary  layers  (WBBL)  in  the  surf  and  swash  zones  and  to 
improve  predictive  capabilities  for  bed  load  and  suspended  sediment  transport  as  a  function  of 
environmental  parameters,  including  wave  heights,  breaker  characteristics,  sediment  properties,  beach 
slope,  bottom  roughness,  local  water  depth,  wave  frequency  spectra,  and  the  presence  of  low  frequency 
circulations  such  as  along  shore  currents  and  undertow. 

OBJECTIVES 

We  are  presently  focused  on  addressing  two  key  aspects  of  sediment  transport: 

1 . )  to  simulate  coupled  two-phase  flow,  particularly  for  bed-load  dominated  flow  regimes,  utilizing 

three-dimensional  hydrodynamic  direct  numerical  simulations  of  the  turbulent,  wave  bottom 
boundary  layer  and  discrete  particle  modeling. 

2. )  to  assess  the  strengths  and  weaknesses  of  the  individual  model  capabilities  and  address  deficiencies 

as  guided  by  the  lab  and  field  data. 

APPROACH 

The  work  involves  theoretical  and  model  development,  numerical  computations,  and  comparison  with 
laboratory  data.  The  primary  experimental  tools  are  three-dimensional  discrete  particle  and  direct 
numerical  simulation  hydrodynamics  models  of  oscillatory,  turbulent  boundary  layers. 
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WORK  COMPLETED 


The  model  being  developed  here  will  challenge  and  refine  existing  parameterizations  for  bedload  and 
suspended  load  transport  rates.  Fundamental  concepts  used  in  describing  the  phenomena  of  sediment 
transport  such  as  the  reference  concentration,  bed  failure  criterion,  and  the  more  recently  introduced 
concept  of  acceleration-induced  transport  will  be  more  accessible  for  study  with  the  new  model.  It  will 
produce  the  high  level  of  detail  necessary  to  refine  our  present  understanding  of  sediment  transport 
processes  and  clarify  new  directions  in  the  measuring  techniques  needed  to  improve  present  predictive 
capabilities. 

Our  approach  integrates  the  turbulent  boundary  layer  hydrodynamics  model  of  Moneris  and  Slinn 
(2004)  (Figure  1)  with  the  discrete  particle  model  of  Drake  and  Calantoni  (2001)  (Figure  2).  The 
models  are  being  integrated  with  two-way  coupling  between  the  fluid  and  solid  phases.  The 
hydrodynamic  model  uses  control-volumes  slightly  larger  than  individual  grain  sizes.  The  fluid  model 
passes  the  three  instantaneous  Cartesian  components  of  momentum  in  a  turbulent  boundary  layer  under 
oscillatory  free-stream  flows  to  the  particle  locations.  The  momentum  exchange  (updated 
approximately  1000  times  per  second)  is  calculated  using  an  empirical  drag  coefficient.  Global 
momentum  is  conserved,  and  as  a  particle  accelerates  the  local  fluid  velocity  decreases,  or  as  the 
particle  decelerates,  the  local  fluid  velocity  increases.  It  is  not  perfect,  but  it  is,  arguably,  a  rational 
step  towards  increased  realism  in  both  the  hydrodynamic  and  discrete  particle  modeling  approaches. 

There  are  several  initial  complexities  with  coupling  the  models  that  we  are  addressing  in  the 
development  stages.  First  is  matching  the  spatial  scales  of  the  two  models.  To  produce  a  reasonable 
turbulence  field  the  hydrodynamics  model  needs  to  include  at  least  5-10  horizontal  eddies  in  each 
direction  of  the  periodic  domain.  Typical  small  eddy  diameters,  in  simulations  over  smooth 
boundaries  and  sand  ripples  (Barr,  Slinn,  Pierro,  Winters,  2004)  are  around  0.5  cm  to  1  cm.  Therefore, 
we  desired  to  have  a  model  domain  on  the  order  of  2-10  cm  to  allow  the  hydrodynamics  to  be  free  to 
develop  in  a  more  natural  manner.  We  settled  on  a  4  cm  horizontal  length  scale  for  the  domain.  To 
accommodate  this  hydrodynamic  constraint,  the  particle  model  was  optimized  to  allow  more  particles 
in  the  domain.  We  can  simulate  approximately  100,000  particles  for  10-20  seconds  of  simulation  time 
in  about  a  day  or  two  of  CPU  time.  The  second  complicating  issue  is  that  as  the  particles  move  across 
the  domain  they  reside  in  multiple  fluid  control  volumes  and  we  need  to  know  the  mass  of  each 
particle  in  each  control  volume  to  a  high  degree  of  accuracy.  In  order  to  appropriately  attribute  the 
exchange  of  momentum  between  the  fluid  and  particulate  phases,  we  needed  to  know  the  portions  of 
each  particle  in  each  fluid  control  volume  as  a  function  of  time.  This  important  detail  has  been  worked 
out  for  a  particle  crossing  through  a  corner  (as  illustrated  in  Figure  3  below)  it  can  occupy  portions  of 
up  to  8  control  volumes  simultaneously. 

There  is  an  additional  complexity.  There  is  a  new  level  of  sophistication  required  in  the  model  because 
the  fluid  does  not  occupy  all  of  the  space.  This  greatly  complicates  the  mathematical  method  required 
for  solution  of  the  fluid  pressure  field.  For  a  uniform  density  flow  field,  the  coefficients  that  appears  in 
front  of  the  Poisson’s  equation  for  pressure  are  constant  (the  density).  However,  for  the  variable 
density  field  (arising  from  portions  of  the  control  volume  being  filled  with  sand  particles),  the 
coefficients  are  spatially  and  temporally  variable.  They  can  be  calculated,  by  integrating  forward  the 
particle  positions,  but  the  computationally  efficient  direct  pressure  solution  method  that  worked  in  the 
old  fluid  model  must  be  replaced  with  a  slower  multi-grid  pressure  solver  (as  used  by  Barr,  Slinn, 
Pierro,  and  Winters,  2004)  or  an  iterative  technique  (as  used  by  Slinn,  Allen,  Holman,  2000).  Iterative 
techniques  introduce  new  mathematical  uncertainty  to  the  model  solution,  because  they  require 
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convergence  criteria  and  tolerance  levels.  This  aspect  of  the  problem  is  the  main  technical  challenge 
remaining.  We  are  seeking  to  obtain  an  efficient,  stable  algorithm.  Details  of  the  total  model 
algorithm  are  given  below. 

RESULTS 


The  model  equations  to  be  solved  for  the  two-phase  flow  are  presented  (e.g.,  Dong  and  Zhang,  1999; 
Hsu  et  ah,  2004).  The  modified  momentum  equation  is 


P 


d  (su) 
dt 


+  (u  •V)su 


-sWp  +  //V2  (eu )  +  peg  -  freaction  +  sFApplied , 


(1) 


and  the  modified  continuity  equation  is  — +  V*(a/)=  0,  (2) 

where  s=  1  -  c  is  the  local  porosity  of  the  fluid,  c  is  the  particle  concentration  at  the  grid  volume,  p  is 
the  fluid  density,  u  is  the  fluid  velocity,  p  is  the  pressure,  li  is  the  dynamic  fluid  viscosity,  f  is  the 
acceleration  due  to  gravity,  and  freaclion  represents  the  coupling  force  which  is  the  component  of  the 
particle  drag  force  acting  back  on  the  fluid  volume  and  FApplied  is  an  applied  external  pressure  gradient 

that  drives  the  wave  bottom  boundary  layer  flow.  The  pressure  field  is  detennined  iteratively  because 
of  the  non-constant  coefficient  £n+l  in  front  of  V//"  at  the  new  (n+1)  time  level  in  the  3ld  order 
Adams-Bashforth  time  stepping  scheme  using  the  projection  method.  It  is  convenient  to  divide 
through  by  sn+>  and  then  take  the  divergence  of  the  momentum  equations  and  add  them  to  obtain. 
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The  continuity  equation  applied  at  the  n  +  1  time  level - =  VCkw  is  used  to  replace  the  second 

dt 

tenn  on  the  R.H.S.  The  last  3  terms  are  unknown  since  un+\  v"+1,  vv" 1 1  are  still  unknown  until  the 
iteration  has  converged.  So  V2 p  can  be  solved  with  iterative  approximations  to  these.  Finally: 
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Note  thatVP 
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0  in  both  clear  water(f  =  l)  and  the  packed  bed  «  0.3)  zones,  and  these 


unknown  terms  may  converge  rapidly  except  across  the  fluid-particle  interface. 
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We  set  convergence  criteria  in  a  standard  fashion,  J^( 

dPn+l 


n+l,m  n+l,m-l 

U;  ~U. 


at  the  top  and  bottom 


dz 


~P 


A tsn+l  L 


)  <  y  .  For  Boundary  Conditions 
( £W')”  -(sw)  H  I  and  we  want  boundary  conditions  w”  1  =  0  , 


so 


dp 


n+ 1 


dz 
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Ate 


/  \* 

n+1  \£w)b  .  The  governing  equation  for  the  motion  of  a  spherical  particle  is  given  by 


(e.g.  Madsen,  1991) 
du 


P,K  =  (A  -  PfsS  +  -pC[,A\u  -  us\(u  -  us)+  pVs 


Du 

Dt 


+  F„  (3) 


where  ps  is  the  sediment  density,  Vs  is  the  particle  volume,  us  is  the  particle  velocity,  C'd  is  the 
modified  drag  coefficient,  A  is  the  projected  area  of  the  sphere  and  F,  is  the  sum  of  inter-particle 


forces.  The  first  term  on  the  right  hand  side  of  the  equation  represents  the  particle  buoyancy.  The 
second  term  represents  the  particle  drag  force  where  the  drag  coefficient  is  calculated  from  a  fit  to  the 

empirical  drag  law  for  a  sphere,  C'd  =  c*  (24  Re]  1  +  4 Res  2  +  0.4),  where  c  represents  a  modification 

1  -5 

based  on  the  local  particle  concentration  (e.g.  Richardson  and  Zaki,  1954),  c  =  (1  -c  — c2)  2. 


The  third  term  represents  the  applied  horizontal  pressure  force,  which  is  the  driving  force  of  fluid  and 
particle  motion  in  the  model,  written  here  in  terms  of  the  free  stream  fluid  acceleration.  A  number  of 
fluid-particle  interaction  forces  (e.g.  added-mass,  Bassett  history,  Magnus  forces)  have  been  ignored 
with  this  initial  fonnulation.  Consider  the  reaction  term,  freaclion ,  in  (1)  is  the  equal  and  opposite  force 
to  the  particle  drag  force.  The  reaction  term  has  dimensions  of  force  per  unit  volume  in  (1).  The 
particle  drag  found  in  (3)  has  dimensions  of  force.  In  general,  the  reaction  term  in  (1)  may  be  written 


as  freaction  =  — — pC'dA\ii  -  us\(ii  -  u\  The  particle  drag  force  is  assumed  to  be  a  body  force 

2  Vs 

acting  through  the  center  of  mass  of  the  particle.  The  fluid  velocity,  u ,  found  in  (1)  is  the  estimated 
velocity  at  the  center  of  the  sphere.  When  the  fluid  velocity  is  constant  over  a  large  volume  compared 
to  the  particle  then  the  estimate  is  trivial;  the  velocity  of  the  fluid  at  the  center  of  the  particle  is 
estimated  to  equal  the  velocity  of  the  fluid  volume  where  the  center  of  the  particle  is  located.  However, 
for  our  problem  a  sediment  particle  could  span  many  grid  volumes  (of  fluid)  in  the  vertical,  while  at 
most  four  volumes  in  the  horizontal,  the  task  of  estimating  the  particle  drag  force  using  the  empirical 
drag  law  for  a  sphere  while  simultaneously  satisfying  Newton’s  Third  Law  becomes  problematic.  Due 
to  the  stretched  vertical  grid  a  sediment  particle  may  occupy  volume  in  over  20  different  fluid  grid 
points  simultaneously.  Detennining  a  reasonable  estimate  of  the  fluid  velocity  at  the  particle’s  center 
for  use  in  the  empirical  drag  law  is  no  longer  trivial.  There  is  NO  clear  choice  that  seems  obviously 
better  than  the  other  possibilities.  Some  methods  will  be  more  computationally  efficient  than  others, 
but  one  could  argue  that  their  use  sacrifices  accuracy  and  more  importantly,  fidelity  to  the  physics. 


After  careful  consideration,  we  have  chosen  a  method  for  computing  the  velocity  at  the  center  of  a 
particle  for  our  baseline  simulations  that  is  somewhere  in  between  the  extremely  difficult  and  trivial.  In 
practice,  there  are  at  least  two  necessary  constraints  that  need  to  be  satisfied  when  determining  how  to 
compute  the  reaction  term.  The  concentration  of  particles  in  fluid  grid  volumes  must  always  be  less 
than  unity.  More  realistically  the  concentration  should  not  be  allowed  to  exceed  about  0.7.  As  a  result 
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of  this  constraint,  there  will  be  some  trade  off  between  the  scale  of  the  turbulence  resolved  and  the 
largest  particle  allowed  in  the  simulation.  For  baseline  simulations  the  horizontal  grid  spacing  will  be 
uniform  in  both  the  i  and  j  directions.  The  length  of  the  side  of  the  horizontal  square  bounding  a  fluid 
grid  volume  needs  to  be  greater  than  or  equal  to  the  diameter  of  the  largest  sediment  particle.  Simply 
written,  the  first  constraint  is  a  >  0  true  for  all  space  and  time  in  the  simulation. 

The  second  constraint  is  that  Newton’s  Third  Law  must  be  strictly  enforced!  When  a  sediment  particle 
is  much  smaller  than  the  fluid  grid  volume  the  easiest  method  for  obeying  Newton’s  Third  Law  is 
clear.  The  fluid  velocity  used  to  compute  the  drag  force  is  just  the  velocity  of  fluid  at  the  grid  point 
where  the  center  of  the  particle  resides.  The  computed  drag  force  is  generated  entirely  from  a  single 
grid  point  and  is  simply  projected  back  onto  that  grid  point  with  equal  magnitude  and  opposite 
direction.  For  our  configuration  the  particle  may  occupy  volume  in  many  grid  points  simultaneously 
and  it  is  not  possible  to  choose  a  single  grid  point  for  fluid-particle  interactions  without  violating  the 
first  constraint.  Implicitly  assumed  is  that  the  total  mass  of  fluid  and  particles  will  be  conserved  for  all 
time  and  space  in  the  simulations. 


z  =  0.27  cm 


Figure  1.  Turbulent  kinetic  energy  dissipation  rates  in  a  simulation  of  a  turbulent  wave  bottom 
boundary  layer  during  a  phase  of  flow  reversal  The  left  panel  shows  a  vertical  plane,  and  the  right 
panel  shows  a  horizontal  plane  located  2. 7  mm  from  the  boundary  during  flow  transition. 
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Figure  2.  Discret particles  during  a  simulation. 


Figure  3.  The  position  of particles  crossing  control  surfaces  must  be  accounted  for  precisely  to 
accurately  determine  the  mass  of  fluid  and  solid  in  each  control  volume.  A  complex  algorithm  has 
been  developed  to  efficiently  determine  the  location  of  the  mass  of  each  particle  in  motion. 


IMPACT/APPLICATIONS 

Our  models  are  the  most  sophisticated,  detailed  models  of  turbulence  and  sediment  transport  ever 
implemented.  This  approach  should  allow  much  more  detailed  understanding  of  the  complex  physics 
of  two-phase  flows.  The  model  results  will  permit  evaluation  of  bulk  transport  formulas. 

RELATED  PROJECTS 

Modeling  projects  for  the  Sand  Ripple  DRI  are  related. 
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PUBLICATIONS 


None  yet. 
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